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ABSTRACT 

Supermassive black hole binaries (SMBHBs) are expected to emit continuous gravitational waves in the 
pulsar timing array (PTA) frequency band (10~^-10~^ Hz). The development of data analysis techniques aimed 
at efficient detection and characterization of these signals is critical to the gravitational wave detection effort. 
In this paper we leverage methods developed for LIGO continuous wave gravitational searches, and explore 
the use of the J^- statistic for such searches in pulsar timing data. Babak & Sesana 2012 have already used this 
approach in the context of PTAs to show that one can resolve multiple SMBHB sources in the sky. Our work 
improves on several aspects of prior continuous wave search methods developed for PTA data analysis. The 
algorithm is implemented fully in the time domain, which naturally deals with the irregular sampling typical 
of PTA data and avoids spectral leakage problems associated with frequency domain methods. We take into 
account the fitting of the timing model, and have generalized our approach to deal with both correlated and 
uncorrelated colored noise sources. We also develop an incoherent detection statistic that maximizes over 
all pulsar dependent contributions to the likelihood. To test the effectiveness and sensitivity of our detection 
statistics, we perform a number of monte-carlo simulations. We produce sensitivity curves for PTAs of various 
configurations, and outline an implementation of a fully functional data analysis pipeline. Finally, we present a 
derivation of the likelihood maximized over the gravitational wave phases at the pulsar locations, which results 
in a vast reduction of the search parameter space. 



1. INTRODUCTION 

In the next few years pulsar timing arrays (PTAs) are ex- 
pected to detect gravitational waves (GWs) in the frequency 
range 10~^-10~^ Hz. Potential sources of GWs in this fre- 
quency range include supermassive black hole binary sys- 
tems (SMBHB s) (Sesana et al. 11 2008]), cosmic (s uper)strings 
( |01mez et aL|[2 010), inflation (St arobinsky||1979|), and a first 
order phase transition at the QCD scale fCapr ini et al.|20 10). 
The community has thus far mostly focused on stochastic 
backgrounds produced by these sources, but they can mani- 
fest themselves in different ways. Cosmic strings and SMB- 
HBs (in highly eccentric orb i ts) can also produce GW bursts 
dPamour & Vilenkin||200T| [Siemens et aT] [20071 [Leblondj 
et al.|[2009| ) in which the duration of the GW signal is much 
less than the observation time. Sufficiently nearby single 
SMBHBs may produce detectable continuous waves wit h pe- 
riods on the order of years ( Wyithe & Loeb 2003; Sesana] 
et al.[[20Q9t jSesana & Vecchio,,2010j . The concept of a 
PTA, an array of accurately timed milli second pulsars , was 
first conceived of ov er two decades ago fRom ani|1989[ Fos- 
ter & Backer] [ 1990] ). Twenty years later three main PTAs 



are m full operation around the world: the North American 
Nanohertz Observatory for Gravitational waves (NANOGrav; 
Jenet et al.| (|200 9)), the Parkes Pulsar Timing Array (PPTA; 
Manchester) (|2008|)), and the European Pulsar Timing Array 
(EPTA; jJanssen et al.[ ( 12008) )). The three PTAs collaborate 



to for m the International Pulsar Timing Array (IPTA; [Hobbs 
|etaLl ([2010)). 

Prior to the establishment of PTAs, Jenet et al. ( 2004) used 
existing pulsar data to rule out the proposed SMBHB system 
3C66B, a possible source of continuous GWs (the mass of the 
proposed system has since been been lowered significantly 
pguchi et al. 2010) so that it not likely to be detectable with 



current PTAs). In this work, the authors looked for the signa- 
ture of a continuous GW in real pulsar data through the use 
of Lomb-Scargle periodograms and sugges ted a metho d for 
directed searches of known sources. Yardl ey et al.| ( [2010] ) also 
relied on the Lomb-Scargle periodogram to detemiine the sen- 
sitivity of the PPTA to continuous GW sources a s a function 
of GW frequency, van Haast eren & Levin[ (12010) developed 
a bayesian framework aimed at the detection of GW mem- 
ory in PTAs; however, the authors mention that the meth- 
ods presented could be used f or continuous GW sources as 
well. jSesana & VecchiojpOlO ) use an Earth-term only signal 
model to perform a study of SMBHB parameters measurable 



^ Center for Gravitation and Cosmology, University of Wisconsin Mil- 
waukee, Milwaukee WI, 5321 1 



with PTAs using a Fisher matrix approach. jCorbin & Cor 
[msh ( 2010) have developed a Bayesian Markov Chain Monte- 
Carlo (MCMC) data analysis algorithm for parameter estima- 
tion of a SMBHB system in which the pulsar term is taken 
into account in the detection scheme, thereby increasing the 
SNR and improving the accura cy of the GW source location 
on the sky. Recently, Lee et al.[ ( |201 1) have developed param- 
eter estimation techniques baseaon vector Ziv-Zakai bounds 
incorporating the pulsar term and have placed limits on the 
minimum detectable amplitude of a continuous GW source. 
In this work, the authors also propose a method of combining 
timing parallax measurements with single- source GW detec- 
tions to improve pulsar distance measurements. 

In the context of searches for continu ous gravitational 
waves from spinning neutron stars in LIGO, |Taranowski et al.j 
( [1998| ) developed the so-called J^-statistic, the logarithm of 
the likelih ood ratio m aximized over some of the signal pa- 
rameters. [Cutler & Schutz (2005) later generalize d the 
statistic to muIH^detector networks. Very recently, Babak &] 
Sesana ( 2012) have used the J^- statistic to show that in PTA 
data multiple SMBHB sources can be resolved in the sky. In 
this paper we build on this work, and improve on a number 
aspects of prior continuous wave search methods developed 
for PTA data analysis. 
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In Section [2] we review the signal model. In Section [3] we 
discuss the ^statistic in the context of PTA data. Umike 
LIGO implementations of the J^-statistic, our algorithm is im- 
plemented fully in the time domain. This naturally deals with 
the irregular sampling of PTA data and avoids the spectral 
leakage problems that arise when frequency domain methods 
are used on such data. We also account for the timing model: 
fitting out pulsar parameters removes signal power at low fre- 
quencies, at frequencies near 1 yr~^ and 2 yr~^ due to sky lo- 
cation, proper motion, and parallax fitting, and for pulsars in 
binaries, at frequencies near the binary orbital frequency. Our 
approach also naturally incorporates colored noise sources, 
both uncorrected and correlated (for the case when the dom- 
inant noise source is a gravitational wave stochastic back- 
ground). We also develop an incoherent detection statistic 
that maximizes over all pulsar dependent contributions to the 
likelihood. To test the effectiveness and sensitivity of our de- 
tection statistics, in Section[4]we perform a number of monte- 
carlo simulations. We produce sensitivity curves for PTAs of 
various configurations, and show that the performance of the 
incoherent statistic is comparable to the coherent J^-statistic. 
We also present an outline of the implementation of a continu- 
ous wave search pipeline. Finally, in Section [5] we summarize 
our results and conclude with a derivation or the likelihood 
maximized over the gravitational wave phases at the pulsar 
locations, which results in a vast reduction of the search pa- 
rameter space. We leave the exploration of this new statistic 
for future work. 

2. THE SIGNAL MODEL 

In this section we will briefly review the form of the resid- 
uals induced by a non- spinning SMBHB in a circular orbit 
and introduce our notation. The GW is defined as a metric 
perturbation to flat space time defined as 



habit, Q) = el,iSl)Kit, Q) + e^,(mx(t, ft), 



(1) 



where ft is the unit vector pointing from the GW source to the 
SSB, h+, hx and (A = +, x) are the polarization amplitudes 
and polarization tensors, respectively. The polarization ten- 
sors can be converted to the Solar System Bary center (SSB) 
by the following transformation. Following |Wahlquist| ( | 1 987| ) 



we write 



where 



m - 
h - 



el^(Q) = rhamb-harib, (2) 

e^f^iCl) = mahb+naMb, (3) 

-(sin cos (j))x - (sin sin (j))y - (cos 0)z-, (4) 

-(sin0)f + (cos^)3), (5) 

-(cos 6 cos (j))x - (cos 9 sin (j))y + (sin 6)z- (6) 



In this coordinate system, = ti/2-5 and ^ = a are the polar 
and azimuthal angles of the source, respectively, where 5 and 
a are declination and right ascension in usual celestial coor- 
dinates. 

We will write our GW induced pulsar timing residuals in 
the following form: 



sit, Cl) = F^(Cl)As+(t)^F'' (Cl)Asx(t), (7) 



where 



AsA(t) = SA(tp)-SA(te), 



(8) 



and te and tp are the times at which the GW passes the Earth 
and pulsar, respectively, and the index A labels polarizations. 
The functions F^(tt) are known as antenna pattern functions 
and are defined by 



F''(n) = 



1 (m • p) - (n • p) 

2 i+n-p 

(m-p)(h-p) 
1 + Cl-p ' 



(9) 
(10) 



where p is the unit vector pointing from the Earth to the pul- 
sar. Also, from geometry we can write 



tp = te-L(l-hQ- p). 



(11) 



Given these definitions, we can write the GW contributions 
to the timing residuals as (jWahlquist|1987l|Corbin & Cornish 
[20TDI ) 



sin[2(<l>(0 - $o)](l + cos^ t) cos It/j 



Sx(t) = 



Dujity/^ 

2 cos[2(<l>(0 - $0)] cos i sin lij; 

Dujity/^ 

+ 2 cos[2(^(0 - $0)] cos i cos lij) 



sin[2(<l>(0-^o)](l+cos^ 0sin2?/; 



where 



and 



1 



32>/5/3 



ujit)=[uj-^''-^-^M"'t 



-5/3 



-3/8 



(12) 



(13) 



(14) 



(15) 



For reasons that will become clear later, we write the residuals 
for pulsar a in the following form 

rocit,(l) = Sait,(l) + nait) 

4 

= Y,[aiiCi,^^,^)A'^it,e,^,uj^)\ (16) 

-hPa(tXi^^^O,i^,0,(l),UJo,La)-^na(t), 

where ( = M^^^D~^, n^it) is the noise in each pulsar and 

Pa = F\Cl)s^(tp)^F''(Cl)sAtpy (17) 

Hereon we will refer to the summation term as the Earth term 
and p as the pulsar term. We write the combination of chirp 
mass and distance to the binary as one parameter because the 
two can not be disentangled unless there is a measurement of 
/, which we do not consider here. It is customary to label 
the parameters ((,L,^o,ijj) and (O,(j),ujo) extrinsic and intrin- 
sic parameters ( Jaranowski et al. 1998 ), respectively. We then 
define the amplitudes and time dependent basis functions 

ai = ( [(1 +cos^ Ocos$oCOs2?/; + 2cos^sin<l>osin2?/^] 

= -C [( 1 + cos^ i) sin $0 cos 2?/^ - 2 cos t cos $0 sin 2?/;] 
^3 = C [( 1 + cos^ l) cos $0 sin 2?/; - 2 cos t sin <l>o cos 2?/^] 
a4 = -([(l-\- cos^ l) sin ^0 sin 2?/^ + 2 cos l cos ^0 cos 2?/^] 
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and 



a2 =F^(fi)w(0-i/3cos(2$(f)) 



(19) 



Throughout this work we assume that the source is slowly 
evolving (i.e. the phase is independent of the chirp mass) and 
uj(t) ^ uq and ^(t) ^ ujot- 

3. THE LIKELIHOOD FUNCTION AND THE J^-STATISTIC 

Here we will introduce our formalism and derive the likeli- 
hood and J^- statistic (the likelihood maximized over extrinsic 
parameters) for PTAs. We will also discuss the statistics of the 
J^-statistic in the presence and absence of a signal and show 
that we obtain the expected behavior for PTA data. 

3.1. Likelihood 

For a pulsar timing array with M pulsars we define the like- 
lihood function of the noise as multivariate gaussian 



P(n) = 



1 



exp --n^E^^n 



where 



(20) 



(21) 



is the vector of the noise time-series for all pulsars, 
E« ■ 



^21 S«,2 



Sim S2, 



M 



Sim 
S2M 



is the multivariate covariance matrix, and 



^n,i = (n/ii/) 
Sij = {ninj)\. 



(22) 



(23) 
(24) 



are the auto-covariance and cross-covariance matrices of the 
pulsar noise, respectively. It is important to note that in the 
case of uncorrected noise, the off-diagonal cross covariance 
matrices vanish. In order to time pulsars, a timing model is 
fitted out of the pulsar times-of-arrival (TOAs) via a weighted 
least squares fitting routine (Hobbs et al. 2006 ). This proce- 
dure c an be expressed via a data-independent linear operator 
R (see pemorest et al.|2012| for details) so that 



n = Rn, 



where 



R = 



^1 



R 



M. 



(25) 



(26) 



is a vector of matrices the fitting operators for each pul- 
sar, and n is the post-fit noise. We can see the effect of this 
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Figure 1. SMBHB waveforms in two different regimes. Each plot shows the 
waveform before (dotted blue) and after fitting (solid green) for a full timing 
model including spin-down, astrometeric and binary parameters. Top Panel: 
The Earth and pulsar terms at the same frequency. Bottom Panel: The Earth 
term and pulsar term at different frequencies. 

fitting procedure on the waveforms in Fig. [T] where the wave- 
form is changed, quite significantly, from its pre-fit form. It is 
straightforward to show that the likelihood for the fitted n is 



where 



S^=(nn^)=R(nn^)R^. 
The fitted residuals can therefore be written as 
r = R(s + n) = s + n, 

where 



(27) 
(28) 
(29) 

(30) 



are the residual data and signal template for each pulsar, re- 
spectively. We can therefore write the likelihood of the data r 
given some signal template s 

p(r|s) = -^=L= exp (^-i(r-g)^S;:Hr-§)) . (31) 

We define the inner product for two time vectors x and y 
using the noise covariance matrix as 

(x|y)=x^S^'y (32) 
In this notation we can write the log of the likelihood ratio as 



In A = In 



P(r\s) 
P(r\0) 



= (r|s)-^(s|s). 



(33) 



It is worth pointing out that finding the inverse of T^n is com- 
putationally intensive. Aside from it being a very large ma- 
trix, the fitting procedure results in loss of degrees of freedom 
in the data which makes T>n singular. Inverting this matrix 
therefore requires singular value decomposition. 

In most realistic scenarios we can assume that the off- 
diagonal cross-covariance matrices are small and expand the 
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inverse of E g. 22 in a Neumann series (see Eq. 72 of Anholm 
|et al.|2QQ9] for oetails). In the simulations shown later in the 
paper we will assume that any correlated noise is much less 
than the uncorrelated part, thus we treat T^fi as a block diago- 
nal matrix of the auto-covariance matrices for each pulsar. 

3.2. The Earth-term -statistic 

We now analytically maximize over the extrinsic parame- 
ters (C, i-, ^0, V^) in the signal model. A very similar calcula- 
tion was first done by Jaran owski et al.|(|l 998) in the context 
of LIGO, subsequently by Cornish & Porter ( 2007 ) in the con- 
text of LISA, and very recently by Babak & Sesana (2012 ) in 
the context of pulsar timing. For clarity, here we review this 
calculation in the notation introduced above. For this calcula- 
tion we treat the pulsar term as a noise source and write our 
signal template in the form 

4 



1=1 



where 



(35) 



Later we will explain the circumstances under which it is safe 
to drop the pulsar term. We can now write the log-likelihood 
as 

In A = a/(r|AO- ^(A'"|AO^i/^?; = ^//N^ - ^M%ay. (36) 

Maximizing the log-likelihood ratio over the four amplitudes 
ai gives 



= = K5\-]-M^aj5\-\M^ai5] 
dak 2 ^ ' 2 ^ 



(37) 



yielding the maximum likelihood estimators for the four am- 
plitudes 

^.• = M,-,-N^ (38) 

where Mij = (M^^)~^ Substituting these back into the likeli- 
hood results in the -statistic 



(39) 



The statistics of 2J> are a with 4 degrees of freedom and a 
non-centrality parameter p^. It is straightforward to show that 
the expectation value is 



(2J-,)=4 + p2 



= 4 + (s|s) + 2(p|s) + (p|A0M,-/p|AO, 



(40) 



where p is the functional form of the pulsar term and the sec- 
ond two terms in are due to the fact that we h ave o nly 
included the Earth term in our templates s. In Figs. |2(c)| and 
2(d)| we can see that the probability distribution functions of 
2 J> follow the expected distributions in the absence and pres- 
ence of a signal. While only the intrinsic parameters are for- 
mally searched over, it is also possible to get estimates of the 



maximize d extrinsic parameters by constructing the following 
quantities ( [Cornish & Porter|2007| ): 



A+ = \/ {cii +a4)^ + (a2-ai,)^ 

+ V (^1 ~^4)^ + (^2+^3)^7 

Ax = ^/ (ai +a4)^ + (a2-ai,)^ 
-^{ai-aA)^^ (a2-\-a3)^ 



and 



A=A+ + JAl+A^ 



It is then possible to recover the maximized parameters 



L = COS 



-Ay 



-tan 
2 

-tan" 

Ac 
T' 



/ A+a4-Ay,Ai \ 
\Axa3+A+a2 J 
-(A^ai -A+a4) 



(41) 
(42) 

(43) 

(44) 
(45) 
(46) 
(47) 



where c = sgn(sm2tlj). It is interesting to examine the case of 
one pulsar. In this case, Eq. [37] has no solution because the 
matrix M is singular. The reason for this is that it is incorrect 
to write the residuals in the form of Eq. [16] with four degrees 
of freedom. For one pulsar, the signal has only two degrees 
of freedom: an amplitude and a phase, or equivalently, two 
unknown amplitudes, thereby making the maximization over 
four independent amplitudes an ill-posed problem. Thus, at 
least two pulsars are needed to solve Eq. It should be 
noted that it is straightforward to generalize tms statistic to N 
GW sources, we will simply have 4N indepe ndent amplitudes 
instead of just 4 (see |Babak & Sesana|2012| for more details). 
However, for simplicity in this work we will deal with just 
one GW source. 

3.2.1. Justification for dropping the pulsar term 

There are two cases in which the pulsar term is truly negli- 
gible to the J^g-statistic and can be dropped from the analysis 
with no change in the statistics. 

The first is the astrophysically likely scenario in which the 
evolution of the GW frequency is such that the Earth and pul- 
sar terms are in different frequency bins (see e.g. Figure 2 of 
,Sesana & Vecchio 2010) . At the frequency of the Earth term 
the signal will up coherently. The pulsar term signals, 
even if they all happen to be at the same frequency, will not 
because they have different phases that depen d on th e the pul- 
sar distances. This effect is illustrated in Fig. 1 3 (a) [ where the 
reference distribution has 4 degrees of freedom and non- 
centrality parameter (s|s). 

The second case is the less astrophysically likely scenario 
in which the Earth and pulsar term lie in the same frequency 
bin. In this case, there is still a phase difference between the 
Earth and pulsar terms. We expect that for a large number 
of pulsars the pulsar term signals will cancel because they all 
have different phases. We can see from Fig. |2(d)| that for a 
moderate number of pulsars (M = 20 in this case) the pulsar 
phases do not completely cancel and our measured values of 
the statistic are higher than expected with just the earth 
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Figure 2. Histograms and expected probability distribution functions of ITp and ITe in the absence and presence of a signal for 20 pulsars. Each simulation 
was done with the search parameters fixed and 1000 realizations of white gaussian noise, (a): distribution of ITp in the absence of a signal, (b): distribution of 
ITp in the presence of a signal with non-centrality parameter . (c): distribution of ITe in the absence of a signal, (d): distribution of ITe in the presence of a 
signal. The dashed (red) and solid (green) curves have the same meaning as in Fig.|3] 

It is indeed possible to include the pulsar term in our anal- 
ysis if we operate in the low frequency (or low chirp mass) 
regime where the frequency evolution of the source is slow 
enough that the frequency at the Earth and the pulsar are es- 
sentially the same so that the signal is a sum of two sinusoids 
of different phases: the pulsar term and the Earth term. To 
understand this more quantitatively, consider the Taylor se- 
ries expansion of the orbital frequency of Eq. [15] evaluated at 
the pulsar time 



term because the last two terms of Eq. 40 do not sum to zero. 
However, in the case of large M (M ^~jO) the pulsar term 
contributions sum approximately to zero, and again we have a 
distribution with 4 degrees of freedom and non-centrality 
parameter (s|s). 

If we happen to detect a signal that falls into the intermedi- 
ate category mentioned above where M < 50 and some or all 
of the pulsar terms are in the same frequency bin as the Earth 
term, then this will create a bias in the recovered sky loca- 
tion but not in our ability to confidently detect the signal (see 
Fig.l of |Ellis et al.| ( |2012a| ). This is because our detection 
criterion is the fals e ala rm probability. As will be discussed 
in detail in Section |3.4[ the false alarm probability only de- 
pends on the probability distribution function when the signal 
is absent, and we can see from Fig. 2(a)] that ITe follows the 
expected distribution, because it is independent of the signal 
properties. 

3.3. The incoherent T -statistic 



uj(tp) = ujo yl 

^ CJo ( 1 + 



256 



-3/8 



te-L(l-\-ft- p) 



(48) 



From this, we can see that uo(tp)^ cjq when 



T-L(l^Q- 



3/8 



(49) 



6 



Ellis et. al. 



0.035 
0.030 
0.025 
0.020 
0.015 
0.010 
0.005 
O.OOQ 




40 60 
J^-Statistic 



0.025 



0.020 



0.015 



0.010 



0.005 



0.000, 



Full Signal ^ 
J^-statistic 




20 40 60 80 100 120 140 160 
J^-Statistic 



(a) 



(b) 



Figure 3. Probability distribution functions for 2J> in the limits that the pulsar term is negligible, (a): probability distribution function in the limit that all pulsar 
terms lie outside the Earth term frequency bin. (b): probability distribution function in the limit of large M for overlapping Earth and pulsar term frequencies. 
The dashed (red) curve is a distribution with a non-centrality parameter assuming that only the Earth term is present in the data. The solid (green) curve is a 
distribution with non-centrality parameter that takes both the Earth and pulsar term into account. 

Maximizing the likelihood ratio over the 2M amplitude pa- 
rameters biciC^ h V^, *^o, (t^a.O, (/)) gives 



where T is the total observation time. If we consider only one 
intrinsic parameter, ujq, then the template for pulsar a is 

2 

So^it, n) = Y^ bio^iC, i, V^, <l>o, , 0, (t))B'^{t, cjo), (50) 



where 



(51) 



is the pulsar dependent phase. We can now write the pulsar 
dependent amplitudes and basis functions as 



( 1 + cos^ 0(^a + 2?/;)(cos <l>o - cos ^q,) 



+ 2 cos b{F^ sin 111) - cos 2?/^)(sin $o - sin (/)c) 



(52) 



( 1 + cos^ i)(F^ cos 2 + F^ sin 2 V^)(sin - sin ) 



- 2 cos i(F^ sin 2?/^ - F^ cos 2?/^)(cos - cos (pa ) 



and 



Bl(t)=—sm(2uJot) 

^0 



(53) 



(54) 



(55) 



where, again, uq is the angular orbital frequency of the 
SMBHB. The log-likelihood ratio is 



In A: 



M 



(56) 



ainA 

dbkp 



=o=E 



(57) 



which yields the solution for the maximum likelihood estima- 
tors of the 2M amplitudes 



(58) 



Putting the amplitude estimators back into the likelihood ratio 
we obtain the J"„-statistic 



=1 



(59) 



It is straightforward to then show that 2Fp follows a dis 
tribution with 2M degrees of freedom and non-centrality pa 
rameter and that 

{2:Fp)=2M+p^ 



= 2M+(s|s) 



(60) 



where = (s|s) is the optimal signal-to-noise ratio (see Fig. 
[2]). Note that this is an incoherent detection statistic since it 
involves sum of the squares of the data, whereas the Earth- 
term J> -statistic is coherent since it involves the square of the 
sum of the data. 

It is worth pointing out that for the case of white gaussian 
noise, the -statistic is the time domain equivalent t o the 
weighted pow er spectral summing technique studied in Ellis 
et al.|(|2012al). For colored gaussian noise the statistic is the 



time domain equivalent to a weighted power spectral sum- 
ming technique with frequency dependent weights. Another 
feature of this detection statistic is that it does not only apply 
to the low-frequency limit. If we work in the high frequency 
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regime where the Earth and the pulsar terms are in different 
frequency bins, we can drop the pulsar term and arrive at the 
exact same maximized likelihood function. In this case the 
pulsar dependence of the amplitudes bia comes from the an- 
tenna pattern functions the not the pulsar phase. However, 
many of the justifications for dropping the pulsar term men- 
tioned in the previous section do not apply in this case since 
the statistic is incoherent. We find that this detection statistic 
will often pick out the pulsar term frequency over the Earth 
term frequency because the residuals of Eq. [16] scale like 
uj(t)~^^^ and the pulsar term will always be at an equal or lower 
frequency than the Earth term frequency due to the geometri- 
cal delay in Eq. [TT] For this system of equations we have 2M 
equations and 6 + M unknowns, so if we have 6 or more pul- 
sars we can solve for the all the parameters (C, ^, V^, ^o,^,^) 
along with the pulsar phases (j)a . 

3.4. False alarm probability and detection statistics 

Here we review the false alarm and detection probability 
distribution functions both when the intrinsic parameters are 
known and unknown. O ur discussion follows closely tha t of 
IJaranowski eTHl i fTWSl ) and |Jaranowski & Kr61akl ( |2n05] ). In 
the case of k now n ex trinsic parameters, we have shown in 



Sections 3.2 and 3.3 that the statistics 2J^ and 2J^p follow 



X distributions with 4 and AM degrees of freedom, respec- 
tively, when the signal is absent. It was also shown that the 
aforementioned statistics follow a non-central with non- 
centrality parameters p and p, respectively, when the signal is 
present. 

Therefore, the probability distribution functions /?o and pi 
when the intrinsic parameters are known and when the signal 
is absent and present, respectively, are 



^n/2-l 



{n/2-l)\ 

(2J^)in/2-l)/2 
1 



(2^)(»/2-i)/2 



: exp 



2 



(61) 



(62) 



where n is the number of degrees of freedom, 4/2- 1 is the 
modified Bessel function of the first kind and order ^/2- 1, 
and ni^ p for J^p and p for The false alarm probability Pp 
is defined as the probability that exceeds a given threshold 
J^Q when no signal is present. In this case, we have 



/^f(^o) = / P^{:F)d:F = exp(-7o) ^ ^ • (63) 



The probability of detection Pd is the probability that J' ex- 
ceeds the threshold To when the signal-to-noise ratio is k,: 



Pd(J^o,i^)= / pi(T,K)dT, 



(64) 



however; we do not deal with the detection probability in this 
work. Our detection criterion is based on the false alarm prob- 
ability. 

We now turn to the more realistic problem of calculating 
the false alarm probability when the intrinsic parameters are 
not known. A detailed derivation and description is given in 
Jaranowski & Krolak ( 2QQQj , here we will simply review the 



result. The probability Pp that exceeds J^o in one or more 
cells is given by 



pJ(Jo)=1-[1-/^f(^o)] 



Nc 



(65) 



where Nc is the number of independent cells in parameter 
space. The number of independent c ells can be calculated 
via geometrical methods described in [Jaranowski & Krol^ 
([2000j) and references therein. 
Here we will make the following approximations. For our 
statistic we will set Np to be equal to the number of inde- 
pendent frequency bins defined by the Nyquist frequency. For 
our Fe statistic, we will set Nc to be equal to the number of 
templates used in the search. In general the number of inde- 
pendent templates and the number of independent cells will be 
quite different. However, since we only have a three dimen- 
sional parameter space and use a nested sampling algorithm 
to conduct the search (thereby reducing the number of tem- 
plates in low likelihood regions of parameter space), setting 
the number of templates equal to the number of independent 
cells is a reasonable assumption. 

4. PIPELINE, SENSITIVITIES, AND IMPLEMENTATION 

In this section we will test the J> and J^p statistics on real- 
istic simulated data sets. First, we will outline our detection 
pipeline, then we will briefly describe our simulated data sets 
and test the ability to confidently detect the signal and recover 
the injected intrinsic parameters. Finally, we perform monte- 
carlo simulations to produce sensitivity curves for PTAs of 
various configurations and sensitivities. 

4.1. Detection Pipeline 

The only inputs to our detection pipeline are the ephemeris 
file (typically called a "par" file) and TOA file (typically 
called a "tim" file) for each pulsar. The steps in the pipeline 
are as follows: 

1. Use the standard pulsar timing package TEMP02 
( Hobbs et al.|2Q06| ) to form the residuals for each pul- 



sar. 

2. Use Tempo 2 plugin to out put the design m atrix for 
each pulsar (see Chapter 15 of |Press et al.|1992| for more 
details ). Then con struct R from the design matrices fol- 
lowing |Demorest| ( |20Q7 ). 

3. Use a maximum lik elihood eigenvalue decomposition 
method described in [Ellis et"aL] ( |2Q12b] ) to make an es- 
timate of H^. Note that the cross terms in Eq. [22] are 
expected to be small, so we will ignore them rar this 
work. 



4. Follow the methods described in Sees. 13.21 and [331 to 
construct the detection statistics and search the relevant 
parameter space. If using the -statistic we simply 
grid up the frequency space for the search. If using the 
J>-statistic we use the nested sampling package, Multi- 
Nest (Feroz et al.|2009| ) to search the three dimensional 
parameter space. 

5. Output the maximum value of the detection statistic 
and number of templates used and compute the relevant 

Here we set our 



false alarm probability using Eq. 65 



false alarm probability threshold to 10~^. If the false 
alarm probability corresponding to our maximum value 
of is greater than 10~^ then we claim a detection. 
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6. Use the maximum likelihood estimators to find the ex- 
trinsic parameters (using Eqs. |44}j47l), and construct 



the posterior probability distribution to find the intrin- 
sic parameters by sampling the maximized likelihood 
(Eq. [39]). As mentioned above, when using the J^p 
statistic, one could use numerical techniques to obtain 
estimates of the extrinsic parameters. 

7. Use the maximum likelihood values of the intrinsic and 
extrinsic parameters to construct Gaussian prior distri- 
butions and carry out parameter estimation on the the 
full 7 dimensional search space, again using MultiNest, 
to get better estimates of SMBHB parameters. 

In this paper we will only conduct steps 1-5 and leave steps 
6 and 7 for future work. Although this work uses simulated 
datasets, nothing in this detection pipeline makes any assump- 
tions about the spacing of the data, or the color of the noise. 

In the absence of a detection we would like to set upper 
limits on the strain amplitude as a function of GW frequency. 
This can be accomplished as follows 

1 . Run the detection pipeline and determine the value of 
the J^- statistic. 

2. For each frequency, choose the value of C correspond- 
ing to a specific strain amplitude. Then inject a 
SMBHB signal with randomly drawn binary orienta- 
tion parameters (cos ^, ?/^, ^o). 

3. Run the detection pipeline again on this injected data 
and measure the value of the J^-statistic. 

4. Keep the value of C fixed and perform a given number 
of injections with different binary orientation parame- 
ters (1000, for example) and determine the fraction of 
J^-statistic values that is larger than the value measured 
in the original data. 

5. Repeat steps 2-4 until the strain amplitude is such that 
95% of the injections give a value of the J^- statistic that 
is larger than the original value. 

6. Record this value and repeat steps 2-5 at each fre- 
quency. 

4.2. Simulated data sets 

For this work we use a simulated pulsar timing array with 
sky locations drawn from uniform distributions in cos^ and (j). 
All pulsars are assumed to have a distance of 1 kpc and a white 
noise rms of 100 ns with equal error bars. The timespan of the 
observations for all pulsars is 5 years with evenly spaced bi- 
monthly TOA measurements. Each set of residuals has been 
created by fitting a full timing model includi ng spin-down, 
astrometric, and binary parameters (see Edwards et al.||2006| 
for details). As a check in some simulations an uncorrelateJ 
red noise process with a power law spectrum P{f) = Af~^ 
is included in the residuals, which has no effect on our re- 
sults. While these simulated data sets do not include uneven 
sampling or extra fitting procedures like jumps or time vary- 
ing DM variations, they do capture the essence of real timing 
residuals in the quadratic fitting of the spin-down parameters 
and the yearly and half yearly sinusoidal trends due to the 
sky location, proper motion and parallax fitting. Very uneven 
sampling is likely to reduce our sensitivity at higher frequen- 
cies and a detailed study of this problem will be presented in 
future work. 



4.3. Implementation of the detection statistics 

Here we will test our detection statistics on mock data sets 
with injected SMBHB GW signals in the presence of white 
and red gaussian noise. We will focus primarily on the 
statistic since, as we will show, it is a more robust detection 
statistic. Then, we will implement a procedure to produce 
an upper limit on the GW strain ampl itude as a function of 
frequency for a simulated NANOGrav |Demorest et al.|20T2 ) 
array and plausible SKA arrays. 

Fig. |4] shows the posterior probability distributions of the 
intrinsic search parameters for simulated SMBHB signals in 
the presence of 100 ns white noise (Fig. |4(a)| ) and uncor- 



related red noise with amplitude A = 4.22 x 10~ s and 
7 = 4.1. The two cases do have different realizations of the 
white noise, however, we can see that the J> -statistic does a 
very good job of determining the frequency and sky location 
of the source. In general, the J>-statistic is more robust than 
the J^p statistic because it produces estimates of the sky loca- 
tion as well as the frequency, which is very important when 
looking for electromagnetic counterparts. 

It is possible to produce a sensitivity curve by a method that 
is similar to what we use to set upper limits. In this case we 
use simulated data with a given level of noise an d no signal 
present. We follow the method presented in Sec. 4.1 except 
we now look for strain amplitude that gives a false alarm prob- 
ability that is higher than our threshold (10~^ in our case) in 
95% of realizations for each frequency. For clarity, we define 
the strain amplitude as 



h = 2 



D 



(66) 



where /gw = ujq/it. This amplitude comes from the overall 
scaling factor that results in differentiating Eq. fT2]and[T3]with 
respect to time. For simplicity and speed we have simplified 
this method for our sensitivity plots. Instead of performing a 
search at each frequency, we simply evaluate the Fe and Fp 
statistics at the values of the injected parameters. The purpose 
of these sensitivity plots is to illustrate the overall features of 
the different detection statistics and to give order of magni- 
tude estimates of expected sensitivity for real data. 

We have produced various sensitivity curves for both the J> 
and Fp statistics in Fig. [5] The three scenarios that we look at 
are a 17 pulsar simulatea NANOGrav array in which we use 
the real sky location and timing models of the NANOGrav 
pulsars, and simulated PTAs with 25 and 100 pulsars at ran- 
dom sky locations. The loss in sensitivity at GW frequencies 
of lyr~^ and 2yr~^ are due to the fitting of the pulsar's sky 
location and proper motion, and parallax, respectively. It is 
important to note that the sensitivity curves for the J> and Fp 
statistics in the 17 and 25 pulsar cases, respectively, are very 
similar. Conversely, for the case of 100 pulsars the Fe statis- 
tic is more sensitive by a factor of ~ 2 for almost all frequen- 
cies. This is due to the different scaling relations of the statis- 
tics vs. the number of pulsars {Fe cx y/M while Fp (xM^^^^). 
However, the plot shows that the -statistic is more sensitive 
at lower frequencies and the J> -statistic is more sensitive at 
higher frequencies. There are two effects that contribute to 
this. The first is a result of our simulation and stems from 
the fact that we assume that for a given frequency, the max- 
imum value of the statistic is at the injected sky location. 
However, for low frequencies where the Earth and pulsar term 
are in the same frequency bin this assumption breaks down as 
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Figure 4. Posterior probability distribution functions for sky location and orbital frequency for a network SNR=14 injection with and without red noise. Here 
we have used a PTA with 25 pulsars. The vertical lines indicated the injected parameters and the contours are the one, two and three sigma contours, (a): 100 
ns white noise, (b): 100 ns white noise and uncorrected red noise with amplitude A = 4.22 x 10~^^ s~^-^^ and 7 = 4.1. We see that the sky location and orbital 
frequency have all been recovered at the one-sigma level in both cases. 



the sky location will be biased (see e.g. [Ellis et aL]|2012a| ). 
The second effect is one inhere nt to our detection statistics 
themselves. As discussed in Sec. 13.31 the J^„-statistic has dif- 



ferent meanings in the low and high frequency regimes. In the 
low frequency regime, it effectively contains the entire signal 
(Earth and pulsar terms), and in the high frequency regime 
it only contains the Earth term piece since the pulsar terms 
are out of that frequency bin. This distinction results in a 
different scaling relation for the ratio of TelTp. In the low 
frequency case the statistic scales coherently but it only 
has approximately half of the signal, whereas, the J^^-statistic 
scales incoherently but has the full signal. Therefore, the ratio 
scales as thus the incoherent method will do better 
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Figure 5. Sensitivity curves for the J> and Tp statistics for different PTA 
configurations (all pulsars have 100 ns residuals). The blue and green lines 
are the sensitivity curves for the J> and Tp statistics, respectively, for a sim- 
ulated NANOGrav PTA. The red and cyan lines are the sensitivity curves for 
the J^e and J^p statistics, respectively, for a simulated PTA with 25 pulsars. 
The magenta and yellow lines are the sensitivity curves for the J> and Tp 
statistics, respectively, for a simulated PTA with 100 pulsars. 



for M < 16. Conversely, in the high frequency regime, both 
statistics contain only half of the signal and the ratio scales as 
M^/^. Therefore, the coherent statistic will do about a factor 
of 2 better than the incoherent method for M > 16. 

5. SUMMARY AND OUTLOOK 



In this work we have adapted the standard J^- statistic (I Jara-| 
|nowskietal.|1998| ) to act as a detection statistic for continuous 
wave searches in realistic PTA data. We have also developed 
an incoherent detection statistic that maximizes over all pul- 
sar contributions to the likelihood. Both of these detection 
statistics are implemented in the time domain to avoid spec- 
tral leakage problems associated with Fourier domain meth- 
ods applied to irregularly sampled data. These methods take 
the pulsar timing model fitting into account and have been 
generalized to account for both correlated and uncorrelated 
colored noise. Most of our analysis relies on dropping the 
pulsar term from our signal model as it will not add coher- 
ently. We have justified the use of this approximation in most 
astrophysically likely scenarios. It was shown that both of 
the detection statistics follow well known distributions in 
the presence and absence of GW signals and therefore have 
well defined false-alarm probabilities. We have shown that 
the statistic can not only confidently detect a GW signal 
but can also determine the sky location and frequency of the 
source to relatively high accuracy in the presence of white 
and colored gaussian noise. A realistic implementation of a 
fully functional continuous GW pipeline starting from basic 
pulsar timing data and methods for computing upper limits 
on the strain amplitude were outlined in detail. Finally, we 
have used simulated data sets of various PTA configurations 
to produce sensitivity curves for our J^-statistics. From these 
sensitivity curves, we have shown that the sensitivity of the 
Te and Tp statistics are very similar for M < 25 pulsars and 
that the statistic becomes more sensitive for M > 25 and 
for higher frequencies. 



As was shown in Ellis et al. ( 2Q12a l), explicitly searching 
over the pulsar distances or somewEat equivalently, the GW 
phases at the pulsar locations (in the low frequency regime), is 
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computationally prohibitive for M > 5. A statistic that could 
maximize over these GW phases would greatly reduce the pa- 
rameter space of the search, while still preserving the SNR 
of the full signal. The implementation of such an algorithm 
will be the subject of future work. However, we will give the 
derivation here. From Eq. [16] in the low-frequency limit we 
can write the signal in the following form 



M 



Sait) = [(cos^a " l)Sij sm^ aSij] A\ (67) 



i=0 



where <I>q, = cjLq,(1 + • p^), d = at and A' are defined in Eqs. 
Island [19] respectively, and the matix 



s - 



"0 -1 0" 

10 

1 

10 



(68) 



After some algebra, the log-likelihood ratio of Eq. |33]can be 
written as 



with 



M 

In A = ^ [Z7(cos^$cK-sin^<l>c) + ccos$a 

q;=1 

+ Jsin^c^ +/sin<l>Q,cos<l>Q,] , 

c = W^ai+M'^aiaj 
d^Nlsija^ 



(69) 



(70) 

(71) 
(72) 
(73) 



/ = -M'^eijaia" 
where Mq, and are defined by the following relations 

^2' = (<l^i) (74) 
< = (r,K). (75) 

Maximizing the log-likelihood with respect to the pulsar 
phases <l>/3, we obtain 



^^-^ =/(cos^<l>/3-sin^^/3) + 2/?cos<l>/3sin<l>/3 
-csin<l>/3+ Jcos<l>/9 = 0. 



(76) 



/3 ■ 



Setting X = cos<l>^, this expression reduces to a quartic equa- 
tion of the form 



= (4f + 16/?^)/ + (4fd + ^cb)x^ 
+ (c^ -4f-l 6b^)x^ + (-2/J - 8c/?)x 



(77) 



which is guaranteed to have at least one unique solution. This 
maximization results in a monumental reduction in the param- 
eter space that needs to be searched. It takes on the order of 
~ 10^^ templates just to cover the pulsar phases (Ellis et al.| 
|2012a). In practice, we could construct the various quantities 
No,, a, b, c, d, and /, solve Eq. 77 numerically to find 
the maximum likelihood estimators for all the pulsar phases. 
Substituting these solutions back into our likelihood Eq. [69] 
still leaves us with the problem of searching over a 7 dimen- 
sional parameter space (since the amplitudes a depend on 4 



parameters (C,^,^o,V^) and the basis functions A depend on 
3 parameters (^,^,a;o)). We note, however, that this can be 
easily handled with a Markov chain Monte-Carlo (MCMC) 
or nested sampling algorithm. 

Looking to the future, the pipeline outlined in this paper 
will be used to analyze real pulsar timing residuals and, in 
the absence of a detection, construct upper limits on the strain 
amplitude as a function of frequency. We will also further 
develop and test our likelihood maximized over the GW phase 
at the pulsar on both simulated and real data. We will also 
begin to generalize the methods discussed in the this paper to 
deal with eccentric signal models. 
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data analysis working group for their comments and support, 
especially Jim Cordes, Paul Demorest, Rick Jenet, Andrea 
Lommen, Delphine Perrodin, Sam Finn, and Joe Romano. We 
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design matrices. This work was partially funded by the NSF 
through CAREER award number 0955929, PIRE award num- 
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